﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature.Numerics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GMix
{
    class IsentropicCompressionSoneFunction : ISoneFunction
    {
        #region Private Fields
        Action<double[], double[]> m_selectedFunction;
        double m_t0;
        double m_p0;
        double m_fixedPressure;
        double m_fixedTemperature;
        double m_fixedCompressionRatio;
        double m_s0;
        double m_expDeltaS;
        CompositionObject m_composition; 
        #endregion

        internal PolytropicGasRatios Configure(CompositionObject composition, Pressure pressure)
        {
            Initialize(composition);
            m_fixedPressure = pressure;
            m_selectedFunction = this.FixedPressureFunction;
            double gamma = composition.CalcGamma(composition.BaseTemperature.Value);
            var ratios = PolytropicGasRatios.CalcIsentropicRatios(composition.BasePressure.Value, pressure, gamma);
            return ratios;
        }

        internal PolytropicGasRatios Configure(CompositionObject composition, CompressionRatio compressionRatio)
        {
            Initialize(composition);
            m_fixedCompressionRatio = compressionRatio;
            m_selectedFunction = this.FixedCompressionRatioFunction;
            double gamma = composition.CalcGamma(composition.BaseTemperature.Value);
            var ratios = PolytropicGasRatios.CalcIsentropicRatios(compressionRatio, gamma);
            return ratios;
        }

        internal PolytropicGasRatios Configure(CompositionObject composition, Temperature temperature)
        {
            Initialize(composition);
            m_fixedTemperature = temperature;
            m_expDeltaS = Math.Exp(
                        (composition.CalcEntropy(m_fixedTemperature) - m_s0) / Constants.Rgas
                    );
            m_selectedFunction = this.FixedTemperatureFunction;
            double gamma = composition.CalcGamma(composition.BaseTemperature.Value);
            var ratios = PolytropicGasRatios.CalcIsentropicRatios(composition.BaseTemperature.Value, temperature, gamma);
            return ratios;
        }

        private void Initialize(CompositionObject composition)
        {
            m_composition = composition;
            m_t0 = composition.BaseTemperature.Value;
            m_p0 = composition.BasePressure.Value;
            m_s0 = composition.CalcEntropy(m_t0);
        }

        #region ISoneFunction Members

        public int Nva
        {
            get { return 3; }
        }

        public void Compute(double[] x, double[] y) { m_selectedFunction(x, y); }

        private void FixedPressureFunction(double[] xarray, double[] y)
        {
            IsentropicCompressionVector x = (IsentropicCompressionVector)xarray;

            double fxp = m_fixedPressure;
            double t = x.Temperature;
            double p = x.Pressure;
            double cr = x.CompressionRatio;
            double s = m_composition.CalcEntropy(t);

            y[0] = s - m_s0 - Constants.Rgas * Math.Log(fxp / m_p0);
            y[1] = cr - fxp / m_p0 * (m_t0 / t);
            y[2] = p - fxp;
        }

        private void FixedCompressionRatioFunction(double[] xarray, double[] y)
        {
            IsentropicCompressionVector x = (IsentropicCompressionVector)xarray;

            double fcr = m_fixedCompressionRatio;
            double t = x.Temperature;
            double p = x.Pressure;
            double cr = x.CompressionRatio;
            double s = m_composition.CalcEntropy(t);

            y[0] = s - m_s0 - Constants.Rgas * Math.Log(t / m_t0 * fcr);
            y[1] = p - m_p0 * t / m_t0 * fcr;
            y[2] = cr - fcr;
        }

        private void FixedTemperatureFunction(double[] xarray, double[] y)
        {
            IsentropicCompressionVector x = (IsentropicCompressionVector)xarray;

            double fxt = m_fixedTemperature;
            double t = x.Temperature;
            double p = x.Pressure;
            double cr = x.CompressionRatio;
            double s = m_expDeltaS;

            y[0] = p - m_p0 * m_expDeltaS;
            y[1] = cr - p / m_p0 * m_t0 / t;
            y[2] = t - fxt;
        }


        #endregion
    }


    class IsentropicCompressionSoneSolver : SoneNewton<IsentropicCompressionSoneFunction>
    {
        internal IsentropicCompressionSoneSolver() : base( new IsentropicCompressionSoneFunction() )
        {
            base.SetAbsoluteConvergenceRadius(
                    IsentropicCompressionVector.Create(50.0, 0.0, 3.0)
                );
            base.SetRelativeConvergenceRadius(
                    IsentropicCompressionVector.Create(0.1, 0.1, 0.5)
                );
            base.SetAbsoluteTolerance(0.0);
            base.SetRelativeTolerance(1.0e-8);
        }
    }
}
